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Highly supercooled liquids with soft-core potentials are studied via molecular dynamics simu- 
lations in two and three dimensions in quiescent and sheared conditions. We may define bonds 
between neighboring particle pairs unambiguously owing to the sharpness of the first peak of the 
pair correlation functions. Upon structural rearrangements, they break collectively in the form of 
clusters whose sizes grow with lowering the temperature T. The bond life time rt,, which depends on 
T and the shear rate 7, is on the order of the usual structural or a relaxation time Tc in weak shear 
7Ta <C 1, while it decreases as I/7 in strong shear jTa 3> 1 due to shear-induced cage breakage. 
Accumulated broken bonds in a time interval (~ 0.05ri,) closely resemble the critical fluctuations of 
Ising spin systems. For example, their structure factor is well fitted to the Ornstein-Zernike form, 
which yields the correlation length ^ representing the maximum size of the clusters composed of 
broken bonds. We also find a dynamical scaling relation, rt, ~ valid for any T and 7 with z = 4 
in two dimensions and 2 = 2 in three dimensions. The viscosity is of order for any T and 7, 
so marked shear-thinning behavior emerges. The shear stress is close to a limiting stress in a wide 
shear region. We also examine motion of tagged particles in shear in three dimensions. The diffusion 
constant is found to be of order Tj^" with u = 0.75 ~ 0.8 for any T and 7, so it is much enhanced in 
strong shear compared with its value at zero shear. This indicates breakdown of the Einstein-Stokes 
relation in accord with experiments. Some possible experiments are also proposed. 
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I. INTRODUCTION 

Particle motions in supercooled liquids are severely re- 
stricted or jammed, thus giving rise to slow structural 
relaxations and highly viscoelastic behavior Re- 
cently much attention has been paid to the mode coupling 
theory, || which is a first analytic scheme describing 
onset of s ow structural relaxations considerably above 
Tg. There, the density fluctuations with wave numbers 
around the first peak position of the structure factor are 
of most importance and no long range correlations are 
predicted. For a long time, however, it has been ex- 
pected that rearrangements of particle configura- 
tions in glassy materials should be cooperative, involving 
many molecules, owing to configuration restrictions. In 
other words, such events occur only in the form of clus- 
ters whose sizes increase at low temperatures. In normal 
liquid states, on the contrary, they are frequent and un- 
correlated among one another in space and time. Such an 
idea was first put forth by Adam and Gibbs jsj, who in- 
vented a frequently used jargon, cooperatively rearranging 
regions (CRR). However, it is difficult to judge whether 
or not such phenomenological models are successful in 
describing real physics and in making quantitative pre- 
dictions. 

Molecular dynamics (MD) simulations can be power- 
ful tools to gain insights into relevant physical processes 
in highly supercooled liquids. Such processes are often 
masked in averaged quantities such as the density time 



correlation functions. As a marked example, we men- 
tion kinetic heterogeneities observed in recent simula- 
tions [§-18|. Using a simple two dimensional fiuid, Mu- 
ranaka and Hiwatari visualized significant large scale 
heterogeneities in particle displacements in a relatively 
short time interval which was supposed to correspond to 
the P relaxation time regime. In liquid states with higher 
temperatures. Hurley and Harrowell | pd| observed simi- 
lar kinetic heterogeneities but the correlation length was 
still on the order of a few particle diameters. The char- 
acterization of these patterns has not been made in these 
papers. Recently our simulations on model fiuid mix- 
tures in two and three dimensions [|T3|-p^ have identified 
weakly bonded or relatively active regions from breakage 
of appropriately defined bonds. Spatial distributions of 
such regions resemble the critical fluctuations in Ising 
spin systems, so the correlation length ^ can be deter- 
mined. It grows up to the system size as T is lowered, 
but no divergence seems to exist at nonzero tempera- 
tures [ |l3| , p^ -pi| . Donati et al. have observed string-like 
clusters whose lengths increase at low temperatures in 
a three dimensional binary mixture [p^ . In addition, 
Monte Carlo simulations of a dense polymer by Ray and 
Binder showed a significant system size dependence of 
the monomer diffusion constant, which indicates hetero- 
geneities over the system size [Q. 

Most previous papers so far are concerned with near- 
equilibrium properties such as relaxation of the density 
time correlation functions or dielectric response. ^From 
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our point of view, these quantities are too restricted or 
indirect, and there remains a rich group of unexplored 
problems in far-from-equilibrium states. For example, 
nonlinear glassy response against electric field, strain, 
etc. constitutes a future problem p^ . In this paper we 
apply a simple shear flow = jy in the x direction and 
realize steady states |^ . The velocity gradient 7 in the y 
direction is called the shear rate or simply shear. We shall 
see that it is a relevant perturbation drastically changing 
the glassy dynamics when 7 exceeds the inverse of the 
structural or a relaxation time r^. As is well known, 
increases dramatically from microscopic to macroscopic 
times in a rather narrow temperature range Gen- 
erally, in near-critical fluids and various complex fluids, 
nonlinear shear regimes are known to emerge when 7 ex- 
ceeds some underlying relaxation rate p3| . However, in 
supercooled liquids, it is unique that even very small 
shear can greatly accelerate the microscopic rearrange- 
ment processes. Similar effects are usually expected in 
systems composed of very large elements such as colloidal 
suspensions. 

Though rheological experiments on glass-forming flu- 
ids have not been abundant, Simmons et al found that 
the viscosity 77(7) = Uxy/j exhibits strong shear-thinning 
behavior. 



^(7) = ??(0)/(l+7^^), 



(1.1) 



in soda-lime-silica glasses in steady states under shear 
p^-p6t, Tjj being a long rheological time. After applica- 
tion of shear, they also observed overshoots of the shear 
stress before approach to steady states. Our previous 
reports have treated nonlinear rheology in su- 

percooled liquids in agreement with these these exper- 
iments. Interestingly, similar jamming dynamics has be- 
gun to be recognized also in rheology of foams 27-29]] 
and granular materials |]30|] composed of large elements. 
Shear-thinning behavior and heterogeneities in configura- 
tion rearrangements are commonly observed also in these 
macroscopic systems. 

As a closely related problem, understanding of 
mechanical properties of amorphous metals such as 
Cu57Zr43 has been of great technological importance 
plt-^5|. They are usually ductile in spite of their high 
strength. At low temperatures T < 0.6 ~ O.TTg local- 
ized bands (< 1/Ltm), where zonal slip occurs, have been 
observed above a yield stress. At relatively high tempera- 
tures T > 0.6 ~ O.TTg, on the other hand, shear deforma- 
tions are induced homogeneously (on macroscopic scales) 
throughout samples, giving rise to viscous flow with 
strong shear thinning behavior. In particular, in their 3D 
simulations Takeuchi et al. |^,^ followed atomic mo- 
tions after application of a small shear strain to observe 
heterogeneities among poorly and closely packed regions, 
which are essentially the same entities we have discussed. 
Our simulations under shear in this paper correspond to 
the homogeneous regime at relatively high temperatures 
in amorphous metals. 



Another interesting issue is as follows. Several ex- 
periments have revealed that the translational diffusion 
constant Z? of a tagged particle in a fragile glassy ma- 
trix becomes increasingly larger than the Einstein- Stokes 
value Des — ksT /2'iTria with lowering T, where 77 is the 
zero-shear) viscosity and a is the diameter of the particle 
^-38|. In particular, the power law behavior D cx rj~'^ 
with V = 0.75 was observed at sufficiently low tempera- 
tures 1^^. Thus D/Des increases from of order 1 up to 
order 10^ ~ 10"^ in supercooling experiments. Further- 
more, smaller probe particles exhibit a more pronounced 
increase of the ratio Drj/T (x D/Des with lowering T 
pTf . It is generally believed that rj is proportional to the 
a relaxation time Tq or the rotational relaxation time 
l/Drot for anisotropic molecules (Drot being the rota- 
tional diffusion constant) [ ^6p^j39| ]. Therefore, individ- 
ual particles are much more mobile at long times t > Ta 
than expected from the Stokes-Einstein formula. In a 
MD simulation on a 3D binary mixture with N = 500 in 
3D the same tendency was apparently seen despite 
of their small system size. Very recently, in a MD simula- 
tion in a 2D binary mixture with N — 1024, Perera and 
Harrowell have observed clear deviation from the linear 
relation D en Ta where Ta is obtained from the decay of 
the time correlation function as in our case in Section 6 
pf . We will examine this problem in a much larger 3D 
system with N — N1 + N2 — 10** generally in the presence 
of shear, where the viscosity and the diffusion constant 
both vary tremendously in strong shear (7 > I/tq). 

The organization of this paper is as follows. In Sec. II 
our model binary mixtures and our simulation method 
will be explained. In Sec. Ill bonds among particle pairs 
will be introduced at distances close to the first peak po- 
sition of the pair correlation functions. Breakage of such 
bonds will then be followed numerically, which exhibits 
heterogeneities enhanced at low temperatures. Their 
analysis will yield the correlation length in Sec. IV. Rhe- 
ology of supercooled liquids will be studied in Sec.V. 
These effects were briefly reported in our previous re- 
ports p"^-p^. In Sec. VI new results on motion of tagged 
particles will be presented. 



II. MODEL AND SIMULATION METHOD 

We performed MD simulations in two dimensions (2D) 
and three dimensions (3D) on binary mixtures composed 
of two different atomic species, 1 and 2, with Ni = N2 = 
5000 particles with the system volume V being fixed. Pa- 
rameters chosen are mostly common in 2D and 3D. They 
interact via the soft-core potential ||9|-|l5|Hl|- 



Vapir) = e{aai3/r)'^^, Gap = ^((Ta + ap), (2.1) 

where r is the distance between two particles and a,l3 = 
1,2. The interaction is truncated at r = 4.5cri in 2D 
and r = 3(71 in 3D. The leapfrog algorithm is used to 
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integrate the differential equations with a time step of 
0.005ro, where 

TO = {m^alltfl\ (2.2) 

The space and time are measured in units of a\ and tq. 
The mass ratio is mijmx = 2, while the size ratio is 

a2/ai = 1.4 (d==2), = 1.2 (d = 3), (2.3) 

where d is the space dimensionality. This size difference 
prevents crystallization and produces amorphous states 
in our systems at low temperatures. 
We fixed the particle density at 

n = 0.8/ct^, (2.4) 

where n — n\ + ni is the total number density. The sys- 
tem linear dimension is i = 118 in 2D and L — 23.2 
in 3D. Then our systems are highly compressed. In 
fact, the volume fraction of the particles may be esti- 
mated as i^{p\n\ + a\n2) = 0.93 in 2D and as |7r(crfni -I- 
(j\n2) = 0.57 in 3D, where overlapped regions are dou- 
bly counted. In such cases, according to the Henderson 
and Leonard theory p3|-^5[| , our binary mixtures may be 
fairly mapped onto one-component fluids with the soft- 
core potential with an effective radius defined by 

a, (3=1, 2 

where x\ = ni/n and X2 = n2ln = 1 — xi are the com- 
positions of the two components and are 1/2 in our case. 
As in the one component case the thermodynamic state 
is characterized by a single parameter (effective density) , 

T,ff=n{elkBTY'''a%f. (2.6) 

For example, Bernu et al. confirmed that the equi- 
librium pressure p may be well fitted to the scaling form, 
p/nksT -19i 6 + 6.848(re//)'', at all xi in 3D. In Ta- 
bles I and II we list T^ff chosen in our simulations to- 
gether with the corresponding scaled temperatures and 
pressures in 2D and 3D. Our pressure data excellently 
agree with the above scaling form for 3D. 

We here introduce the pair correlation hmctions (Jq^ (r) 
by 

{ha{r)h0{O)) = Uanpgapir) + naSapS{r), (2.7) 
where 

""W -I]'5(r-r,j) (a = 1,2) (2.8) 

i 

are the number densities in terms of the particle positions 
raj (a = 1, 2, j = 1, • ■ • , N/2). The time dependence is 
suppressed for simplicity. In a highly compressed state 
the inter-particle distances between the a and (3 particles 
are characterized by 



ea0=(Tai3{e/kBTy^^^. (2.9) 

The last factor {e/kBTy^^^ represents the degree of ex- 
pulsion or penetration from or into the soft-core regions 
(r < ctq/s) on particle encounters, though it is not far 
from 1 in our case. The one-fluid approximation may be 
justified if the pair correlation functions satisfy 

g„0(r) = G(r/C0,re//). (2.10) 

Namely, gapir) are independent of a, (3, and xi once the 
distance is scaled by £a0 ■ The pressure is then expressed 
as II 

= 6V,r,f f r ds^G{s,r,fj), (2.11) 

Jo ^ 

where (r) — dvap (r) / dr and Vd is the volume of a unit 
sphere, so it is 47r/3 in 3D and tt in 2D. We shall see that 
(2.10) excellently holds around the first peak of the pair 
correlation functions in our simulations. This fairly sup- 
ports the one-fluid approximation, because the soft-core 
potential and the pair correlation functions decrease very 
abruptly for r > £ap and for r < £a0, respectively, and 
the dominant contribution arises from r ~ iafS ~ ctq/j- 

In our systems the structural relaxation time becomes 
very long at low temperatures. Therefore, the annealing 
time was taken to be at least 10^ in 2D and 10"*^ in 3D. 
No appreciable aging effect was detected in the course of 
taking data in various quantities such as the pressure or 
the density time correlation function except for the low- 
est temperature cases, Tg// = 1.4 in 2D and Fg// — 1.55 
in 3D. A small aging effect remained in the density time 
correlation function in these exceptional cases, however. 



TABLE I. Simulations in 2D. 





1.0 1.1 1.2 


1.3 


1.4 


ksT/e 


2.54 1.43 0.85 


0.526 


0.337 


p/nksT - 1 


15.1 22.6 33.5 


50.2 


75.1 


TABLE 11. Simulations in 3D. 




1.15 1.3 1.4 1.45 


1.5 


1.55 


ksT/e 


0.772 0.473 0.352 0.306 


0.267 


0.234 


p/nksT - 1 


18.9 26.7 33.4 37.2 


41.4 


46.3 
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erformed in the absence and 
In the unsheared case 



Our simulations were 
presence of shear flow 47 
(7 = 0) we performed simulations under the microcanon- 
ical (constant energy) condition. However, in the sheared 
case (7 > 0), we kept the temperature at a constant us- 
ing the Gaussian constraint thermostat to eliminate the 
viscous heating effect. No difference was detected be- 
tween the profile-based and profile-unbased thermostats 
p7t , so results with the profile-based thermostat will be 
presented in this paper. Our method of applying shear 
is as follows: The system was at rest for f < for a 
very long equilibration time and was then sheared for 
t > 0. Here we added the average velocity 7?/ to the ve- 
locities of all the particles in the x direction at t = and 
afterwards maintained the shear flow by using the Lee- 
Edwards boundary condition . Then steady states 
were realized after a transient time. In our case shear 
flow serves to destroy glassy structures and produces no 
long range structure. 



III. PAIR CORRELATIONS AND BOND 
BREAKAGE 

A. Pair correlations 



(3.5) 



where Speff = Peff — (Peff)- They are linear combina- 
tions of the usual structure factors, 



(3.6) 



from the definitions (3.1) and (3.2). The temperatures 
in Fig. 3 are common to those in Fig. 2. Note that the 
dimensionless wavenumber q is measured in units of cr^^. 
The Spp{q) has a pronounced peak at g ~ 6 and becomes 
very small (~ 0.01) at smaller q both in 2D and 3D. In 
this sense our systems are highly incompressible at long 
wavelengths. On the other hand, Sxxiq) has no peak 
and is roughly a constant over a very wide q region, sug- 
gesting no enhancement of the composition fluctuations 
and no tendency of phase separation at least in our sim- 
ulation times. 

^From Fig. 3 we may estimate the magnitude of the 
isothermal compressibility Ktx = {dn/dp)Tx/'n. In 
equilibrium it is expressed in terms of the fluctuation 
variances as 



Because of the convenience of visualization in 2D, we 
first present a snapshot of particles at Fe// — lA in 2D 
in Fig.l, which gives an intuitive picture of the particle 
configurations. We can see that each particle is touch- 
ing mostly 6 particles and infrequently 5 particles at dis- 
tances close to (Japi— 1.095~^^Q,/3). Similar jammed par- 
ticle configurations can also be found in 3D, where the 
coordination number of other particles around each par- 
ticle is about 12. Then it is natural that the pair cor- 
relation functions gafsir) {a, f3 = 1,2) have a very sharp 
peak at r = CTq.^, as displayed in Fig. 2 for Fe// = 1.4 in 
2D and T^ff = 1.55 in 3D. Furthermore, the heights of 
these peaks are all close to 7 in 2D and 4 in 3D. This 
confirms the scaling form (2.10) around the first peak. 

We newly introduce a density variable representing the 
degree of particle packing by 

(3.1) 



in terms of which the local volume fraction of the soft- 
core regions is 7r/5e//(r) in 2D and by (47r/3)pe//('") in 
3D. We also consider the local composition fluctuation. 
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dX{r) = -[x2ni{r) - Xin2(r) 



(3.2) 



where xi = X2 — in our case. In Fig. 3 we show the 
corresponding, dimensionless structure factors, 



Sppiq) 



dre^^-''{5Mf{r)5Mfm, (3.3) 



Spx{q) - 'y^" / dre^'l-'-{5M}{r)5X{0)). (3.4) 



ksTKTx = n lim 

g-»0 



Sii{q)S22{.q) ~ Si2{qf 



Sxx{q) 



Sppiq) - Spx{qf / Sxx{q) 



(3.7) 



The first line was the expression in Ref. Q , and the sec- 
ond line follows if use is made of (3.1) and (3.2). The 
dimensionless combination nksTRTx is equal to 0.0028 
in 2D and 0.0067 in 3D. If we assume that the adiabatic 
compressibility Ksx — {dn/dp)sx /n is of the same order 
as Ktxi the sound speed c turns out to be of order 10 in 
units of (Ti/ro. 

Our structure factors were obtained by time averag- 
ing over very long times, which are 10^ for 2D and 10^ 
for 3D. However, irregular shapes of Sxx{q) persisted at 
long wavelengths q < 1. Such large scale composition 
fluctuations have very long life times (s> Tq) and are vir- 
tually frozen throughout the simulation. Therefore, we 
admit the possibility that our supercooled states at low 
temperatures might phase-separate to form crystalline re- 
gions on much longer time scales. On the contrary, the 
long wavelength fluctuations of Pe// have much shorter 
time scales, probably they vary on acoustic time scales 
~ 1/cq. 

As is well known, the temperature dependence of the 
static pair correlation functions is much milder than that 
of the dynamical quantities. Similarly, their shear depen- 
dence is also mild even for 7Tq, ^ 1 as long as 7 ^ 1. 
In particular, their spatially anisotropic part is at most a 
few percents of their isotropic part around the first peak 
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positions r = CTq^ in our case. This is consistent with the 
fact that the attained shear stress in our simulations are 
at most a few percents of the particularly high pressure p 
of our systems. Note that the average shear stress axy in 
sheared steady states may be related to the steady state 
pair correlation functions gai3{f) as ||45| ] 

c'xy ^ -^^n^np I drv'^^{r)^^gap{r), (3.8) 

where and ry are the x and y components of the vector 
r connecting particle pairs. The dominant contribution 
here arises from the anisotropy at r = (Jai3- 

B. Bond breakage 

Because of the sharpness of the first peak of gapif) 
in our systems, we can unambiguously define bonds be- 
tween particle pairs at distances close to aafi in the ab- 
sence and presence of shear. Such bonds will be broken 
on the structural (a) relaxation time, because the bond 
breakage takes place on local configurational rearrange- 
ments. We define the bonds as follows. For each atomic 
configuration given at time a pair of particles i and j 
is considered to be bonded if 

r^j{ta) = |rj(to) - rj(io)| < Moap, (3.9) 

where i and j belong to the species a and (3, respectively. 
We have set Ai = 1.1 for 2D and 1.5 for 3D. The resul- 
tant bond numbers between a and [3 pairs, Nta/s, are 
related to the first peak structure of gapir) as follows. 
We consider the coordination number 1^0.(3 of P particles 
around a a particle within the distance Aiaap p^ , 

i^af} =n[3 drgapir) ~ Cnpa'^p, (3.10) 

where C is about 5 in 2D and 12 in 3D. Then we simply 
have 

Nbaa = ^NaVaa {o. = 1, 2), Nbi2 = ]^NiVi2 + ^N2l^21- 

(3.11) 

In 2D at T^ff = 1.4, we find i^u = 2.19, = V21 = 2.54, 
and V22 = 3.41, which are consistent with the bond num- 
bers, Nbii = 5514, iVbii = 13135, and Nb22 = 84 36, 
counted in a simulation. In 3D at T^ff = 1.55, these 
numbers are — 5.57, 1^12 — V21 — 6.90, 1^22 ~ 8.30, 
which are again consistent with Nbu — 13925, Nbn — 
34476, and Nb22 = 20 7 44 in a simulation. We stress 
that our bond definition is insensitive to Ai, owing to 
the sharpness of the first peak, as long as it is somewhat 
larger than 1 and smaller than the second peak distances 
divided by aap- 

After a lapse of time At, a pair is regarded to have 
been broken if 



rij{tQ + At) > A2craf3 (3.12) 

with A2 — 1.6 for 2D and 1.5 for 3D. This definition of 
bond breakage is also insensitive to A2 as long as A2 > Ai 
and A2(7af3 is shorter than the second peak position of 
(r) ■ We have followed the relaxation of the total sur- 
viving (unbroken) bonds Nbond{At) from the initial num- 
ber 

Nbond{Q) = Nbll + Nbl2 + Nb22 (3.13) 

to zero with increasing At. No significant difference has 
been found between the bond breakage processes of the 
three kinds of bonds, 1-1, 1-2, and 2-2, so we consider 
their sum only. We define the bond breakage time by 

Nbondin) = NbondiO)/e. (3.14) 

The relaxation is not simply exponential at low temper- 
atures, apparently because of large scale heterogeneities 
composed of relatively weakly and strongly bonded re- 
gions. If we fit Nbondi^t) to the stretched exponential 
form, Nbond{At) ~ exp[— (At/rh)° ], the exponent a' is 
close to 1 at relatively high temperatures but is consid- 
erably smaller than 1 at the lowest temperatures (for 
example, a' ~ 0.6 at T^ff = 1.55 in 3D). 

In Fig. 4 we show the bond breakage time Tb = Tb{T) in 
the absence of shear as a function of the temperature. It 
grows strongly with decreasing the temperature. As will 
be shown in (6.8) in Sec. Ill, the bond breakage time Tb 
is proportional to the a relaxation time obtained from 
the decay of the self part of the time correlation function 
Fs{q,t) at g = 27r. The shear dependence of the bond 
breakage time Tb = Tb{'j) is also of great interest. As 
shown in Fig. 5, the bond breakage rate 1/Tb{'^) consists 
of the thermal breakage rate l/rfc(0) strongly dependent 
on T and a shear-induced breakage rate proportional to 
7. It is expressed in the simplest conceivable form, 

l/Tbi^)^l/niO) + Ab^, (3.15) 

where Ab = 0.57 in 2D and 0.80 in 3D. In the strong 
shear condition ^Tb{0) > 1, jump motions are induced 
by shear on the time scale of I/7. We shall see that the 
bond breakage occurs more homogeneously with increas- 
ing shear. Therefore, it is natural that, when the strain 
7 = jAt reaches 1, a large fraction of bonds have been 
broken by shear. 



IV. HETEROGENEITY IN BOND BREAKAGE 

Following the bond breakage process we can visualize 
the kinetic heterogeneity without ambiguity and quanti- 
tatively characterize the heterogeneous patterns. In Fig. 6 
we show spatial distributions of broken bonds in a time 
interval of [to, to + 0.05rfc] in 2D, where about 5% of the 
initial bonds defined at t = to have been broken. The 
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dots are the center positions R.y = ^{ri{to) + rj{to)) of 
the broken pairs at the initial time Iq. The broken bonds 
are seen to form clusters with various sizes. The hetero- 
geneity is marked in the glassy case (b) with T^/f = 1.4 
and 7 = 0, whereas it is much weaker for the liquid case 
(a) with Teff — 1 and 7 = 0. The bond breakage time Tb 
is 17 in (a) and 5x10"* in (b). In (c) we set 7 = 0.25 xlO"^ 
and Fe// = 1.4 with Tb = 32 ~ I/7. The heterogeneity 
is known to become much suppressed by shear, while its 
spatial anisotropy remains small. Notice that even in nor- 
mal liquids bond breakage events frequently occur in the 
form of strings involving a few to several particles, obvi- 
ously because of the high density of our system. In glassy 
states such strings become longer and aggregate forming 
large scale clusters. In 3D we also observe string-like 
jump motions in accord with Ref. JlTf and aggregation 
of such strings at low temperatures. 

In Fig. 7 we write the broken bonds in two consecutive 
time intervals, [^07*0 + O-OSt^] and [ig + O.OSt;,, Iq + O.lrb] 
at Teff = 1.4 and 7 = 0. The clusters of broken bonds 
in the two time intervals mostly overlap or are adjacent 
to one another. This demonstrates that weakly bonded 
regions or collectively rearranging regions (CRR) follow 
complex space-time evolution on the scales of ^ and Tf,. 
We do not know its evolution laws but will encounter a 
dynamical scaling law between ^ and tj, in (4.4) below. 

We define the structure factor Sb{q) of the broken 
bonds as 



We also notice that Sb{q) is insensitive to the temper- 
ature at large q, so from the OZ form (4.2) we find 



Sbio) - e 



(4.3) 



The clusters of the broken bonds are thus very analogous 
to the the critical fluctuations in Ising spin systems. In 
fact, small scale heterogeneities with sizes £ in the region 
1 £ <C ^ are insensitive to the temperature. The re- 
lation (4.3) is analogous to the relation, x in 
Ising spin systems between the magnetic susceptibility 
X = lim^^o S{q) and the correlation length ^ near the 
critical point. Here S{q) is the spin structure factor and 
77 is the Fisher critical exponent (^ 1 in 3D). 

Obviously, ^ represents the order of the maximum 
length of the clusters. However, Adam and Gibbs ^ 
intuitively expected that the minimum size of CRR in- 
creases as exp(const./(T — To)) on lowering T towards 
Tq. It has also been discussed as to whether or not 
there is an underlying thermodynamic phase transition 
at a nonzero temperature Tq in highly supercooled liquids 
|]l9|-^ . ^From our data we cannot detect any divergence 
of ^ at a nonzero temperature, although this is not con- 
clusive due to the finite size effect arising from ^ ~ L. 

Furthermore, as in critical dynamics, we have con- 
firmed a dynamical scaling relation between the bond 
breakage time Tb and the correlation length ^, 



Sbiq) 



1 



exp(iq • Hi 



(4.1) 



where the summation is over the broken pairs, Nb is 
the total number of the broken bonds in a time inter- 
val [to, to -l- At], and the angular average over the direc- 
tion of the wave vector has been taken. Furthermore, we 
have averaged over 5-50 Sb{q) data calculated from se- 
quential configurations of broken bonds. Fig. 8 displays 
the resultant Sb{q) after these averaging procedures on 
logarithmic scales at several Fe// without shear. The 
enhancement of Sb{q) at small q arises from large scale 
kinetic heterogeneities growing with increasing Fg/ / both 
in 2D and 3D. /^From a plot of \/Sb{q) versus q^ in our 
previous reports we already found that Sb{q) can be 
nicely fitted to the Ornstcin-Zernikc (OZ) form : 



Sb{q) = ^6(0)/(l 



(4.2) 



The correlation length ^ is determined from this expres- 
sion. It grows up to the system length at the lowest 
temperatures and is insensitive to the width of the time 
interval At as long as it is considerably shorter than the 
bond breakage time Tb Q. The agreement of our Sb{q) 
with the OZ form becomes more evident in the plots of 
Sb{q) / Sb{0) versus q^ in Fig. 9, in which all the data col- 
lapse onto a single OZ master curve both in 2D and 3D. In 
particular, in 3D the deviations are very small, although 
^ L for low T and small 7 in our case. 



n = Ae, 



(4.4) 



where z = 4 in 2D Q and z = 2 in 3D. The coefficient A 
is independent of T^ff and 7 chosen in our simulations, 
as shown in Fig. 10 (a) and (b). Notice that the data 
points at the largest <^ in Fig. 10 are those at zero shear 
for each ^eff- At present we cannot explain the origin 
of these simple numbers for z. We may only argue that 
z should be larger in 2D than in 3D because of stronger 
configurational restrictions in 2D. It is surprising that 
(4.4) holds even in strong shear 7'T{,(0) ^ 1, where the 
correlation length is independent of T and is determined 
by shear as 



-l/z 



(4.5) 



In Fig. 10b for 3D, however, we notice ^ > L at Teff = 
1.50 and 1.55 for weak shear. At present we cannot assess 
infiuence of this finite size effect. 

In a zeroth order approximation, therefore, the kinetic 
heterogeneities are characterized by a single parameter, 
^ or Tb, owing to the small space anisotropy induced by 
shear in our systems. The shear rate 7 is apparently 
playing a role similar to a magnetic field h in Ising spin 
systems. Thus, 7 and T are two relevant external pa- 
rameters in supercooled liquids, while h and the reduced 
temperature (T — Tc)/Tc are two relevant scaling fields 
in Ising systems. 
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V. SUPERCOOLED LIQUID RHEOLOGY 

We next examine nonlinear rheology in our fluid mix- 
tures in supercooled amorphous states. We first display 
in Fig. 11 the shear-dependent viscosity rj^j) (in units of 
ETo/crf) versus 7 in steady states at various Teff in 2D 
and 3D. This rhcological behavior is similar to those in 
the experiments p4|-p6[ . The viscosity is much enhanced 
at large Te// (low T) and at low shear, but it tends to 
be independent of T at very high shear. Remarkably, 
glassy states exhibit large non- Newtonian behavior even 
when 7 is much smaller than the microscopic frequency 
1/ro = 1, whereas such large effects are expected to ap- 
pear only for^~ I/tq in normal liquids far from the 
critical point 

In Fig. 12 we demonstrate that the viscosity 77(7) — 
(^ap/i is determined solely by the bond breakage time 
Tf,(7) in (3.15) as 



(5.1) 



where and 77s are 0.34 and 6.25 in 2D, and 0.24 and 
2.2 in 3D, respectively. Because the linearity 77 oc Tf, 
is systematically violated at small Tf,, the presence of 
the background viscosity rjB independent of T^ff and 
7 may be concluded. Note that the effective exponent 
{j/ri){dri/dj) remains about —0.8 in Fig. 11. As well as 
the kinetic heterogeneities, steady state rheology is deter- 
mined only by a single parameter, Tb or ^. This suggests 
that a sheared steady state can be fairly mapped onto a 
quiescent state with a higher temperature but with the 
same ^. 

Substitution of (3.15) then yields 



^7(7) - A 



(5.2) 



We will argue to derive the above behavior intuitively. 
Supercooled liquids behave as solids against infinitesimal 
strain on time scales shorter than rf,(0) even if the tem- 
perature is considerably above the so-called glass tran- 
sition temperature. Fluid-like behavior is realized only 
after the bond breakage processes. It is natural that the 
viscosity is of order Tf,(0) in the linear regime. This is 
usually justified from the time correlation function ex- 
pression for the viscosity in terms of the xy component 
of the stress tensor In strong shear, on the other 

hand, the bond breakage occurs on the time scale of I/7. 
Upon each bond breakage induced by shear, the particles 
involved release a potential energy whose maximum is 
e. There should be a distribution of e^, but let us assume 
er ^ £ for simplicity. It is then instantaneously changed 
into energies of random motions (and probably sounds) 
supported by the surrounding particles. The heat trans- 
port is rapid in this dissipative process. Because of this 
and also because of the background thermal motions su- 
perposed, we have not detected clear temperature inho- 
mogeneities like hot spots around broken bonds in our 
simulations. The heat production rate is estimated as 



Q ^ ne/Tb{-y) ~ ne'y, 



(5.5) 



where n is the number density. Because Q is related to 
the viscosity by Q = cTxyj — '7(7)7^, we obtain 



77(7)7 ~ ne, 



(5.6) 



in high shear, so crum ^ ne. Due to disordered particle 
configurations, however, it is natural to consider a distri- 
bution of the released energy e^, which will explain the 
viscosity behavior at lower shear. Such a distribution 
was calculated for a model foam system in shear flow by 
Durian E^. 



This form coincides with the empirical law (1.1) by 
Simons et al. p^ , p5[ . Fig. 13 shows that the ratio 
(^7(7) ~ '7s)/('7(0) — ris) can be fitted to the universal 
curve 1/(1 + Abx) with x — juiO) independently of 
Fe// both in 2D and 3D. In strong shear 7T;,(0) ^ 1, 
we have the temperature- independent behavior 77(7) ^ 
{A,-i/Ai,)/'j + rjB, which is evidently seen in Fig. 11. If 
the background viscosity is negligible, a constant limit- 
ing stress follows as 



Ar,/Ab 



which holds for 



l/Tb{0) < 7 < amin/VB - O.I/tq. 



(5.3) 



(5.4) 



Here aiim is 0.59 in 2D and 0.30 in 3D in units of e/af 
and is typically a few percents of the pressure in our 
systems. The upper bound in (5.4) is very large in usual 
glass-forming liquids but should be attainable in colloidal 
systems, while the lower bound can be very small with 
lowering T. 



VI. MOTION OF TAGGED PARTICLES 

In this section we will follow the motion of tagged par- 
ticles in a glassy matrix both in the absence and presence 
of shear in 3D. We will present results only in three di- 
mensions. We first plot in Fig. 14 the self part of the 
density time correlation function for various Fe/ / in the 
usual zero shear condition. 



Fsiq,t) 



1 



^exp[iq • Arj(t)] J, 



(6.1) 



where q — 27r, Arj{t) — rj{t)—rj{0), and the summation 
is taken over all the particles of the species 1. This func- 
tion is proportional to the (incoherent) scattering am- 
plitude from labelled particles. As is well known, this 
function has a plateau at low temperatures (Fe// > 1.45 
in our case), during which the particle is trapped in a 
cage. After a long time the cage eventually breaks, re- 
sulting in diffusion with a very small diffusion constant 



7 



D. In this paper we define the a relaxation time such 
that Fs{q,Ta) = a,t q = 2it. 

We generalize the time correlation function (6.1) in the 
presence of shear flow by introducing a new displacement 
vector of the j-th particle as 

Ar,(t) - r,{t) - 7 / dt'y,{t')e, ~ r,(0), (6.2) 
Jo 

where e^; is the unit vector in the x (flow) direction. 
In this displacement, the contribution from convective 
transport by the average flow has been subtracted, which 
can be known from the time derivative, 

§i^r,{t)^v,it)~jy,it)e.,. (6.3) 

To get clear understanding of the meaning of this sub- 
traction, let us consider a Brownian particle placed in 
shear flow as a simple example. On time scales longer 
than the relaxation time of its velocity, its position r{t) 
obeys 

^^r{t)^jy{t)e, + f{t), (6.4) 

where f{t) is the Gaussian random force characterized 
by iUmAt')) = ^DS^Jit - t') = x,y,z). Then 

the modified displacement vector reads 

Ar{t) = r{t) - 7 / dt'yit')e^ - r(0) = / dt'f (f). 
Jo Jo 

(6.5) 

Here the convective effect does not appear explicitly and 
the diffusion behavior follows as 

(Ar(t)2) = 6Dt. (6.6) 

On the other hand, in the incoherent scattering am- 
plitude, Arj{t) in (6.1) should be taken as the net dis- 
placement rj{t) — rj{0) even in shear flow. li Qx ^0, it 
strongly depends on the thickness of the scattering region 
in the y (velocity gradient jdirection due to a position- 
dependent Doppler ofTcct |g3|,|o|. Only for q^; = 0, it is 
proportional to Fg {q, t) in the above definition. 

Fig. 15 shows Fs{q,t) a,t q = 2tt for various 7 with a 
fixed temperature, ^eff = 1-5 or ksT/e = 0.267 in 3D. 
Comparison of this figure with Fig. 14 suggests that ap- 
plying shear is equivalent to raising the temperature. Re- 
call that we have made the same statement in analyzing 
the bond structure factor S'b(q) and the nonlinear rhcol- 
ogy. Also we may define the shear dependent a relaxation 
time Ta =7-^,(7) by 

F,{q,T^) = e-\ (6.7) 

In Fig. 16 we recognize that Tq, is proportional to the bond 
life time r;, as 



r„ O.lTb. (6.8) 

This relation holds for any Fe// and 7 in our 3D sim- 
ulations. The decay of Fs(q,t) is not exponential for 
large Ta ■ If it is fitted to the stretched exponential form 
exp[— (t/rQ,)"] around t ^ t^, the exponent a is increased 
from values about 0.8 to 1 with increasing 7 as well as 
with raising T. Furthermore, the time correlation func- 
tion (6.1) has turned out to be almost independent of the 
direction of the wave vector q. 

Next it is convenient to analyze the mean square dis- 
placement of tagged particles of the species 1, 

1 

((Ar(i))^) = — ^((Ar,(i))^) (6.9) 

Fig. 17 shows the transition from the ballistic behav- 
ior ((Ar(t))2) ^ 2,{kBT/mi)t^ to the diffusion behavior 
((Ar(t))2) ^ QDt in shear flow at T^ff = 1.55. The ar- 
rows in the figure indicate the a relaxation time Ta(7). 
The crossover occurs around t ^ Ta- Fig. 18 demonstrates 
the surprising isotropy of the statistical distribution of 
Ari(i), where the mean square displacements of the x, y, 
and z components of the vector Ar ^ (t) are separately dis- 
played. We can thus determine D from the mean square 
displacement in addition to in shear flow. Note that 
the X component in Fig. 18 is not the usual mean square 
displacement due to the second term in (6.2). In the 
appendix we will consider the variances of the net dis- 
placement vector rj(t) — Tj(0). 

Fig. 19 shows the shear rate dependences of the vis- 
cosity 77 (~ Ta) and the inverse diffusion constant 
from the linear (7^ lO^'^) to the non-Newtonian regime 
at Teff — 1.55 in 3D, where D is measured in units 
of crf/TQ and 77 in units of eTo/af. We deduce the re- 
lation D^^ ^ 7^'' with ly = 0.75 ^ 0.80 in agreement 
of the experiment ||3^ , which is appreciably milder than 
the viscosity decrease rj ^ Ta ^ 7^^. In Fig. 20 we plot 
D versus rj/ksT (in units of T^/af) obtained for various 
^eff and 7. The Einstein- Stokes formula, which holds 
excellently in normal liquids, appears to be violated in 
supercooled liquids as the other simulations have sug- 
gested ^,0. It is widely believed that this breakdown 
is a natural consequence of the dynamic heterogeneity 
in glassy states p6|-^. Detailed numerical analysis will 
appear in a forthcoming paper. 

In our case rj/ksT changes over 4 decades until ^ 
reaches the system dimension L, whereas it has been 
changed over 12 decades in the experiments p6| , ^ . 
Though the same tendency indicating the breakdown of 
the Einstein-Stokes relation has been obtained in our sim- 
ulation, we should admit that our system size in 3D is 
not yet sufficiently large and our data at Fg// = 1.5 and 
1.55 might be somewhat affected by the system size ef- 
fect. It is worth noting that the Monte Carlo simulation 
of a dense polymer by Ray and Binder shows that 
the monomer diffusion constant decreases with increasing 
the system size. 
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VII. SUMMARY AND DISCUSSIONS 



Most of our findings in this work have been obtained 
from numerical analysis only without first principle 
derivations. Nevertheless, we believe that they pose 
new problems and suggest new experiments. We make 
some discussions mentioning possible experiments below. 

1) Introducing the concept of bond breakage, we have 
succeeded to quantitatively analyze the kinetic hetero- 
geneities in simple model systems, which have been wit- 
nessed by a number of the authors. As shown in Fig. 6, 
strings composed of broken bonds are very frequent and 
they aggregate at low temperatures to form clusters. 
The bond breakage time Tb is related to the correlation 
length ^ as (4.4). In future work we should clarify the 
relationship of our patterns in the a relaxation and those 
by Muranaka and Hiwatari on a much shorter time 
scale. 

2) The weakly bonded regions identified by the bond 
breakage are purely dynamical objects. Large scale het- 
erogeneities have not been clearly detected in snapshots 
of the usual physical quantities such as the densities, the 
stress tensor, the kinetic energy (=temperature), etc. On 
the other hand, in granular matters in shear flow 
stress heterogeneities have been observed optically by 
using birefringent materials. We admit the possibility 
that such stress heterogeneities also exist in supercooled 
liquids but are masked by the thermal fluctuations. Wc 
will check this point in future. 

3) It is of great interest how the kinetic heterogeneities, 
which satisfy the dynamic scaling (4.4), evolve in space 
and time and why they look so similar to the critical fluc- 
tuations in Ising systems in the mean field level. In our 
steady state problem T and 7 are two relevant scaling 
fields, the critical point being located at T = 7 = 0. No 
divergence has been detected at a nonzero temperature 
in our simulations. 



4) In his experiments Fischer |51| has reported large ex- 
cess light scattering with a correlation length ^ (20 ^ 200 
nm) which increases on approaching the glass transition 
from a liquid state. This indicates the presence of very 
large scale density heterogeneities in supercooled liquids, 
which is often called Fischer's clusters. Motivated by 
this effect, Weber et al performed Monte Carlo sim- 
ulations on a dense polymer and found that short range 
nematic orientational order can give rise to enhancement 
of long range density fluctuations. They expected that 
such anisotropic interactions could be the origin of Fis- 
cher's clusters. This suggests that Fischer's clusters do 
not exist in liquids composed of structureless particles. 

5) We have examined nonlinear rheology in glassy states. 
The rheological relations obtained are simplest among 
those consistent with the experiments [pij-psf. The 
mechanism of the non-Newtonian behavior in super- 
cooled liquids is conceptually new and should be further 
examined in experiments such as in colloidal systems in 
glassy states. In particular, polymers should exhibit pro- 



nounced non-Newtonian behavior, as the glass transition 
is approached, even without entanglement. Rheology of 
chain systems remains totally unexplored near the glass 
transition. 

6) In our systems small anisotropic changes of the pair 
correlation functions gaf^ir) near the first peak (~ ctq,^) 
can give rise to the limiting shear stress aum, which is 
3 ^ 5% of the pressure in our case. Note that our systems 
are highly compressed with high pressure. However, the 
pressure needs not be very high in supercooled liquids 
in the presence of an attractive part of the potential. 
Even in such cases, we expect that dum is a few % of the 
shear modulus. This is suggested by the previous work 
on amorphous alloys pl|-|35| , where the yield stress (Jy in 
the inhomogeneous case (in which shear bands appear) 
is known to be 2 ~ 3% of the shear modulus. 

7) Stillinger expected that in fragile glass-forming liq- 
uids shear flow occurs by tear and repair of slipping walls 
separating strongly bonded regions 0. We have not 
observed such localization of slips or jump at least in 
our temperature range. But there might be a tendency 
that broken bonds form surfaces at low temperatures in 
3D, though not conspicuous, which should be checked in 
future. 

8) There is no tendency of phase separation for the pa- 
rameters used. However, there are many cases in which 
the composition fluctuations are enhanced towards the 
glass transition temperature. It is of great interest how 
the two transitions influence each other . It is also 
known that shear flow can induce composition fluctua- 
tion enhancement in asymmetric viscoelastic mixtures, 
when emergence of less viscous regions can reduce the 
effective viscosity We expect that this effect can 
come into play also in supercooled liquids, for example, 
for large enough size ratios or in the presence of small 
attraction between the two components. Experiments 
to detect this effect seem to be promising in colloidal 
systems. 

9) We have introduced the time correlation function 
Fs{q, t) in shear and found its simple relaxation behavior 
in Fig. 15. It coincides with the usual time correlation 
function for = or when the scattering vector is 
perpendicular to the flow direction. Dynamic scatter- 
ing experiments in shear flow would be very informa- 
tive to detect the shear-induced diffusion A direct 
diffusion measurement in sheared supercooled fluids is 
also very interesting, which we will be analyzed in the 
appendix. Though our system size is still too small, 
we have detected a tendency of the breakdown of the 
Einstein- Stokes relation in 3D to obtain D ^ rj^'^ with 

= 0.75 ~ 0.8. 

10) In strong shear the structural relaxation is charac- 
terized by Ta ~ O.lTb ^ 0.1/7 as (6.8). This nonlinear 
effect could be measured as drastic reduction of the rota- 
tional relaxation time by dielectric response or by more 
sophisicated techniques [^,0 from sheared supercooled 
liquids. The same effect is expected for periodic shear 
flow. 
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11) Understanding of transient mechanical response in 
terms of the kinetic heterogeneities is of great impor- 
tance. For example, we have found a stress overshoot 
after application of shear strain in accord with the ex- 
periments 1^,^. We should also understand glassy 
behavior of the complex shear modulus against small pe- 
riodic shear |^^. On these topics we will report shortly. 

12) In our systems we have not yet found essential dif- 
ferences between 2D and 3D except for the difference 
in the value of the dynamic exponent z in (4.4). We 
believe that a large part of essential ingredients of glassy 
dynamics can be understood even in two dimensions. 

13) In a forthcoming paper we will focus our attention on 
jump motions of particles over distances longer than <j\. 
They will be shown to occur heterogeneously in space 
and determine the diffuison constant. These heterogene- 
ity structures are essentially the same as those in the 
bond breakage processes studied in this paper. 



The cross correlation exists between the x and y compo- 
nents as 

((x,(i) - x,{my,(t) ~ ym)) = 7 / dt^G(t^). (A.4) 

Jo 

In the diffusion time regime t >Tawe may set G{t) — 2Dt 
to obtain 

{{x,{t) - xM - jtyof) = 2Dt{l + ^jh'), (A.5) 

{{x,{t) - x,iO)){y,{t) - y,{0))) - D^t'. (A.6) 

Note that D is strongly dependent on 7 in strong shear 
as shown in Fig. 19. Measurements of the above variances 
are very informative. 
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APPENDIX 

Let us calculate the variances among the a;, y, and z 
components, Xj (t) — Xj (0) , yj (t) — yj (0) , and zj (t) — Zj (0) , 
of the net displacement vector rj{t) — rj{0) of the j th 
particle in shear flow. We fix its initial position rj (0) at 
''o = {xo,yo, zq). The average displacement arises from 
the convection as 



(A.l) 



Assuming the isotropy of the subtracted displacement 
(6.2), which is suggested by Fig. 18, we may write the 
variances of the y and z components are 

Git) = mt) - y,m') = {{z.it) - z,(0))2). (A.2) 

The variance of the x component then becomes 



{{xjit)-x,{0)-jtyor) = G{t) + 2Y I dh{t-ti)G{ti). 

'0 

(A.3) 



[1] J. Jackie, Rep. Prog. Phys., 49 171 (1986). 

[2] M.D. Ediger, C.A. Angell and S.R. Nagel, J. Phys. Chem. 

100, 13200 (1996). 
[3] U. Bengtzelius, W. Gotze and A. Sjolander, J. Phys. C 

17, 5915 (1984). 
[4] E. Leutheusser: Phys. Rev. A 29, 2765 (1984). 
[5] G. Adam and J.H. Gibbs, J. Chem. Phys. 43, 139 (1965). 
[6] M.H. Cohen and G.S. Grest, Phys. Rev. B 20, 1077 

(1979). 

[7] F.H. Stillinger, J. Chem. Phys. 89, 6461 (1988). 

[8] K.L. Ngai, R.W. Rendell and D.J. Plazek, J. Chem. 

Phys. 94, 3018 (1991). 
[9] T. Muranaka and Y. Hiwatari, Phys. Rev. E 51, R2735 

(1995). 

[10] T. Muranaka and Y. Hiwatari, Supplement to Prog. 

Theor. Phys. 126, 403 (1997). 
[11] M.M. Hurley and P. Harrowell, Phys. Rev. E 52, 1694 

(1995) . 

[12] D.N. Perera and P. Harrowell, Phys. Rev. E 54, 1652 

(1996) ; J. Non-Cryst. Solids, the proceedings for the 
Vigo meeting (Vigo, June 30- July 11, 1997); 

[13] R. Yamamoto and A. Onuki, J. Phys. Soc. Jpn., 66 2545 

(1997) . 

[14] R. Yamamoto and A. Onuki, Europhys. Lett. 40, 61 
(1997). 

[15] A. Onuki and R. Yamamoto, J. Non-Cryst. Solids, the 
proceedings for the Vigo meeting (Vigo, June 30- July 11, 
1997); the proceedings for the Towa conference (Fukuoka, 
Nov. 4-7, 1997). 

[16] W. Kob, C. Donati, S.J. Plimton, PH. Poole, and S.C. 
Glotzer, Phys. Rev. Lett. 79, 2827 (1997). 

[17] C. Donati, J.F. Douglas, W. Kob, S.J. Phmton, P.H. 
Poole, and S.C. Glotzer, preprint. 

[18] P. Ray and K. Binder, Europhys. Lett. 27 53 (1994). 

[19] C. Dasgupta, A.V. Indrani, S. Ramaswamy and M.K. 
Phani, Europhys. Lett. 15, 307 (1991). 



10 



[20] R.M. Ernst, S.R. Nagel and G.S. Grest, Phys. Rev. B 43, 

8070 (1991). 

[21] S.S. Ghosh and C. Dasgupta, Phys. Rev. Lett. 77, 1310 

(1996) . 

[22] D.J. Salvino, S. Rogge, B. Tigner and D.D. Osheroff, 
Phys. Rev. Lett. 73, 268 (1994) ; S. Roggc, D. Natclson 
and D.D. Osheroff, J. Low Temp. Phys. 106, 767 (1997). 

[23] A. Onuki, J. Phys. C 9, 6119 (1997). 

[24] J.H. Simmons, R.K. Mohr and C.J. Montrose, J. Appl. 
Phys. 53, 4075 (1982). 

[25] J.H. Simmons, R. Ochoa, K.D. Simmons and J.J. Mills, 
J. Non-Cryst. Solids 105, 313 (1988). 

[26] Y. Yue and R. Briickner, J. Non-Cryst. Solids, 180 66 

(1994) . 

[27] T. Okuzono and K. Kawasaki, Phys. Rev. E, 51, 1246 

(1995) . 

[28] D.J. Durian, Phys. Rev. E, 55, 1739 (1997). 

[29] S.A. Langer and A.J. Liu, J. Phys. Chem. B 101, 8667 

(1997) . 

[30] B. Miller, C.O'Hern and R.P. Behringer, Phys. Rev. Lett. 
77, 3110 (1996) ; R. Khosropour, J. Zirinsky, H.K. Pak 
and R.P. Behringer, Phys. Rev. E, 56, 4467 (1997). 

[31] H.S. Chen and M. Goldstein, J. Appl. Phys. 43, 1642 
(1971). 

[32] F. Spaepen, Acta Metall. 25, 407 (1977). 

[33] A.S. Aragon, Acta Metall. 27,47 (1979). 

[34] K. Maeda and S. Takeuchi, Phys. Stat. Sol. 49, 685 
(1978) ; S. Kobayashi, K. K. Maeda and S. Takeuchi, 
Acta Metall. 28 , 1641 (1980); K. Maeda and S. Takeuchi, 
Phil. Mag. A 44, 643 (1981). 

[35] S. Takeuchi and K. Maeda, in Metallic and Semiconduct- 
ing Glasses, ed. A. Bhatnagar, (Trans. Tech. Pub., 1987) 
and references quoted therein. 

[36] F. Fujara, B. Geil, H. Silescu and G. Fleischer, Z. Phys. 
B 88, 195 (1992); I. Chang, F. Fujara, B. Geil, G. 
Heuberger, T. Mangel and H. Silescu, J. Non-Cryst. 
Solids 172-174, 248 (1994). 

[37] M.T. Cicerone, F.R. Blackburn and M.D. Ediger, Macro- 
molecules 28, 8224 (1995); M.T. Cicerone and M.D. Edi- 
ger, J. Chem. Phys. 104, 7210 (1996). 

[38] F.H. StilUnger and A. Hodgdon, Phys. Rev. E, 50, 2064 
(1994). 

[39] D. Richter, R. Frick and B. Farago, Phys. Rev. Lett. 61, 
2465 (1988). 

[40] D. Thirumalai and R.D. Mountain, Phys. Rev. E 47, 
479 (1993). 

[41] D. Perera and P. Harrowell, private communication. They 

have also reproduced the dynamical scaling exponent 

z = 4 in a 2D simulation using a different definition of ^. 
[42] J. Matsui, T. Odagaki and Y. Hiwatari, Phys. Rev. Lett. 

73, 2452 (1994). 
[43] B. Berrm, Y. Hiwatari and J. P. Hansen, J. Phys. C 18, 

L371 (1985); B. Bcrnu, J. P. Hansen, Y. Hiwatari and G. 

Pastore Phys. Rev. A 36, 4891 (1987). 
[44] D. Henderson and P.J. Lenard, Physical Chemistry, an 

Advanced Treatise, ed. D. Henderson (Academic, 1971) 

p. 414. 

[45] J. P. Hansen and I.R. McDonald, Theory of Simple Liq- 
uids (Academic, 1986). 
[46] M.P. Allen and D.J. Tildesley, Computer Simulation of 



Liquids, (Clarendon, Oxford, 1987). 
[47] D.J. Evans and G.P. Morriss, Statistical Meclianics of 

Nonequilibrium Liquids (Academic, New York, 1990). 
[48] J.G. Kirkwood and F.P. Buff, J. Chem. Phys. 19, 774 

(1951). 

[49] H.J.M. Hanley, D.J. Evans, and S. Hess, J. Chem. Phys., 

78 1440 (1983). 
[50] A. Ormki and K. Kawasaki, Phys. Lett. A 72, 233 (1979); 

B.J. Ackerson and N.A. Clark, J. Physique 42, 929 

(1981). 

[51] E.W. Fischer, Physica A 201, 183 (1993). 

[52] H. Weber, W. Paul, W. Kob and K. Binder, Phys. Rev. 

Lett. 78, 2136 (1997). 
[53] A. Sappelt and J. Jackie, Physica A, 240, 453 (1997). 
[54] G. Meier, D. Vlassopoulos and G. Fytas Europhys. Lett. 

30, 325 (1995). 

[55] N. Menon, S.R. Nagel and D.C. Venerus Phys. Rev. Lett., 
73, 963 (1994). 



11 






FIG. 1. A typical particle configuration and the bonds de- 
fined at a given time at T^ff = 1.4 in 2D. The diameters of 

the circles here are equal to aa- The areal fraction of the 
soft-core regions is 93%. A 1/16 of the total system is shown. 
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FIG. 2. The pair correlation functions gatiij) in quiescent 
states as functions of r/aap at Fe// = 1.4 in 2D (a) and at 
Fe// = 1.55 in 3D (b). 
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FIG. 3. The structure factors S{q) defined in (3.3)- (3.5) in 
quiescent states at Fe// = 1.4 in 2D (a) and at Teff = 1.55 
in 3D (b). The dimcnsionlcss wavcnunibcr q is measured in 
units of erf ^. The solid, dashed, and dotted lines correspond 
to Peff—Peff, X — X, and peff—X correlations, respectively. 



FIG. 4. Temperature dependence of the bond breakage 
time Tb{0) at zero shear (•) in 2D and (o) in 3D. The e is 

the potential parameter in the soft-core potentials (2.1). The 
time is measured in units of to in (2.2), so T(,(0) is dimension- 
less. 
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FIG. 5. The normalized bond breakage time T(,(7)/ri,(0) 
versus 7x5(0) for various Fe// in 2D (a) and 3D (b). All the 
data collapse on the curve 1/(1 + Abx) with x = (0). 



FIG. 6. Snapshots of the broken bonds in 2D without 
shear. The system length is llScri. Here Fe// = 1 with 
weak heterogeneity (a), and Fg// = 1.4 with enhanced het- 
erogeneity (b). For 7 = 2.5 x lO"'^ (c), the heterogeneity is 
much suppressed. The flow is in the upward direction and 
the velocity gradient is in the horizontal direction from left to 
right. The arrows indicate the correlation length ^ obtained 
from (4.2). 



14 




FIG. 7. Broken bond distributions in two consecutive time 
intervals, [io, to + O.OSrt] (□) and [to + 0.05t6, to + O.lTb] (•), 
at Teff = 1.4 in 2D. The arrow indicates ^. 
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FIG. 8. Sb{q) versus q on logarithmic scales for various 
Teff at 7 = in 2D (a) and 3D (b). Its long wavelength limit 
is of order as (4.3). 
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FIG. 9. Sb{q)/Sb{0) on logarithmic scales for various Fe// 
and 7 in 2D (a) and 3D (b). The solid line is the Orn- 
stein-Zernike form 1/(1 + x'^) with x = q^. 



FIG. 10. Universal relation between the correlation length 
^(7) in units of cti and the bond breakage time Tb{'j) in units 
of To in (2.2). In 2D (a), the line of the slope 4 is a viewing 
guide and L is the system length. The corresponding 3D plot 
is shown in (b) with the slope being 2. 
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FIG. 11. The viscosity rjl'y) in units of ero/af vesus the 
shear rate 7 in units of I/tq at various Fe// in 2D (a) and 3D 
(b). The data tend to become independent of T^ff at high 
shear. 



FIG. 12. 77(7) vesus s r6(7) for various Fe// in 2D (a) and 
3D (b). The 77(7) is determined by Tb{'y) only irrespectively 

of Te//. 
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FIG. 13. (r;(7) - VB)/iv{0) - Vb) vs 7^(0) in 2D (a) and 
3D (b). The solid curve is 1/(1 + Atx) with x = 7rf,(0). 



FIG. 15. The time correlation function Fs{q,t) at q — 2n 
defined by (6.1) and (6.2) in shear flow, where 7 = 0, 10"'', 
10~^, 10~^, 10~' from right. The temperature is fixed at 
ktT/e = 0.267 (Te// = 1.5). Increasing 7 is equivalent to 
raising T. 
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FIG. 16. The linear relationship between Ta and Tb for var- 
ious Teff and 7 in 3D. The Tb is determined from the bond 
breakage (3.14), and Ta from the decay of the time correlation 
function (6.7). 



FIG. 18. The mean square displacements of the x, y, and 
z components. They are very close to one another even even 
in strong shear 7T(,(0) ^ 1. This demonstrates surprising 
isotropy of the distribution of the displacement vector (6.2). 




FIG. 17. The mean square displacement in sheared states 
at Fe// = 1.5. The shear rate 7 is 0, 10"*, 10"-\ 10"^ 10"^ 
from right. Increasing 7 is equivalent to raising T. The ar- 
rows indicate Ta for each 7. The diffusion (linear) behavior is 
attained at t ~ Ta. 



FIG. 19. Shear rate dependences of the inverse diffusion 
constant D and the viscosity r] at the lowest temperature, 
Teff = 1.55. The slope of is noticeably smaller than 
that of T). 
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ri/kBT 

FIG. 20. The diffusion constant D versus the viscosity rj 
divided by ksT for various Fe// and 7. Here D is measured 
in units of OxTq^ and rj/kBT in units of uf^ro. The sohd 
hne represents the Einstein-Stokes formula D = UbT /2-Rr]a\, 
which well agrees with the numerical data for rj/kBT <10. 
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